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ABSTRACT 


This  thesis  focuses  on  the  numerical  modeling  of  the 
oceanic  double-diffusive  convection,  an  important  small- 
scale  mixing  process  which  is  driven  by  the  two  orders  of 
magnitude  difference  in  dif fusivities  of  heat  and  salt  in 
seawater.  This  study  explores  the  diffusive  regime  of 
double-diffusion  which  is  realized  when  cold  and  fresh  water 
overlies  the  warm  and  salty  water. 

The  aim  of  the  research  was  to  quantify  the  double- 
diffusive  transport  in  both  smooth  gradients  and 
thermohaline  staircases,  and  to  develop  clear  insight  into 
the  origin  of  the  staircases  and  specify  conditions  for 
their  formation.  Based  on  the  numerical  process  modeling,  it 
was  determined  that  the  evolutionary  pattern  of  staircases 
is  controlled  by  the  merging  events  in  which  weak  interfaces 
gradually  erode  and  ultimately  disappear.  To  illustrate 
dynamics  of  these  events,  a  theoretical  framework  -  merging 
theorem  -  has  been  developed.  It  was  numerically  confirmed 
that  the  merging  theorem  predicts  the  time  scale  of  merging 
events  within  the  order  of  magnitude.  The  validity  of  the 
lab  derived  4/3  flux  law  was  tested  and  it  was  determined 
that  its  form  is  consistent  with  the  numerical  results,  but 
the  amplitude  requires  adjustment. 

The  computed  fluxes  from  numerical  experiments  were 
comparable  to  the  diffusive  fluxes  that  were  inferred  from 
the  Beaufort  Gyre  observations (Wilson,  2007)  and  an  order  of 
magnitude  greater  than  the  fluxes  from  earlier  laboratory- 
based  experiments. 

Although  a  large  discrepancy  in  the  values  of  diffusive 
fluxes  exists  in  the  earlier  laboratory  and  field  studies. 


v 


the  present  analysis  suggests  that  the  diffusive  fluxes 
could  play  an  important  factor  in  the  Arctic  heat  budget; 
hence,  future  study  in  this  field  is  recommended. 

This  study  contributes  to  a  better  understanding  of 
global  climate  change,  which  presents  a  new  challenge  to 
national  security.  The  Navy  has  to  be  concerned  with  the 
impact  of  climate  change  on  naval  operations,  specifically 
in  the  Arctic  where  the  melting  polar  ice  cap  may  soon 
provide  a  gateway  across  the  north,  resulting  in  the  opening 
of  shipping  lanes  and  borders. 
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I. 


INTRODUCTION 


A.  DESCRIPTION  OF  DOUBLE -DIFFUSIVE  PROCESS 

1.  Importance  of  Double-diffusive  Process 

The  significance  of  double-diffusive  mixing  and  its 
effect  on  various  large  oceanic  features  has  become  a  major 
driving  force  for  much  of  oceanographic  research  in  this 
field  over  the  last  few  decades.  In  this  thesis  two 
important  aspects  are  explored:  (1)  the  mechanics  of 
diffusive  staircases  and  (2)  the  possibility  that  vertical 
heat  fluxes,  which  are  driven  by  the  diffusive  convection, 
play  a  significant  role  in  the  maintenance  of  the  Arctic 
thermocline  and  melting  rates  of  the  polar  ice  cap. 

2.  Double-diffusive  Process  -  Introduction 

Double-diffusive  convection  may  occur  in  a  stratified 
fluid  if  its  density  is  determined  by  two  or  more 
constituents  with  different  molecular  dif fusivities .  This 
process  is  commonly  observed  in  the  ocean  where  the 
diffusivity  of  heat  is  about  100  times  greater  than  the 
diffusivity  of  salt  (typical  value  of  the  molecular 
diffusivity  of  heat  ~  1.4xl0~7  m2s_1,  and  that  of  salt  ~ 

l.lxlO-9  m2s_1)  .  Since  heat  diffuses  much  faster  than  salt, 
the  resulting  instability  is  vigorous  and  often  leads  to 
elevated  mixing.  Double-diffusive  convection  was  first 
discovered  by  Henry  Stommel  (Stommel  et  al.,1956)  and  Melvin 
Stern  (Stern,  1960). 
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3. 


Salt-finger  Regime  of  Double-diffusive  Process 


Two  modes  of  the  double-diffusive  process  are  known: 
the  salt-finger  regime  and  the  diffusive  convection.  To 
understand  the  salt-finger  mode,  consider  a  case  in  the 
ocean  where  a  warm  and  salty  layer  of  water  overlies  a 
colder  and  fresher  layer  of  water,  the  bottom  layer  having  a 
higher  density  than  the  top  layer.  In  the  absence  of 
molecular  diffusion,  if  a  parcel  of  water  is  displaced 
downward  by  some  force  from  the  lower  density  layer  into  the 
higher  density  layer,  it  tends  to  return  to  its  initial 
position  due  to  buoyancy  force.  However  the  situation 
becomes  very  different  when  diffusion  is  taken  into  account. 
Instead,  due  to  the  higher  diffusivity  of  heat,  the  parcel 
will  become  denser  (and  heavier)  than  the  surrounding  water 
and  will  continue  to  sink  further.  As  a  result,  narrow, 
elongated  vertical  structures  are  formed,  which  are  commonly 
referred  to  as  salt  fingers. 

The  faster  diffusing  component  (heat)  is  stabilizing, 
and  the  slower  diffusing  component  (salt)  is  destabilizing. 
The  regions  favorable  for  the  salt  finger  regime  include  the 
tropics  and  subtropics  where  both  the  temperature  and 
salinity  decrease  with  depth.  According  to  Schmitt  (1990), 
90%  of  the  upper  kilometer  of  the  Atlantic  at  24  degrees 
North  is  favorable  to  salt  finger  formation. 

Other  favorable  regions  are  semi-enclosed  seas  which 
outflow  into  the  ocean,  such  as  the  Mediterranean  Sea  and 
its  outflow  through  the  Straits  of  Gibraltar.  Due  to  the 
exceedance  of  evaporation  over  precipitation  in  the 
Mediterranean,  the  warm  and  salty  waters  of  the 
Mediterranean  flow  over  the  colder  and  fresher  waters  of  the 


2 


Atlantic  and  spread  out  to  a  depth  of  about  1  kilometer.  A 
numerical  example  is  presented  in  Figure  1. 


S 
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Figure  1.  Numerical  simulation  of  salt  fingers.  Red  color 
corresponds  to  high  values  of  salinity;  low  values  are 

shown  in  blue. 

4 .  Diffusive  Convection  Regime  of  Double-diffusive 
Process 

The  second  regime,  the  diffusive  convection,  occurs 
when  a  cold,  fresh  layer  overlies  a  warmer,  saltier  layer. 
Again,  in  order  to  preserve  stability,  the  bottom  layer  has 
greater  density  than  the  top  layer.  Due  to  differing 
dif fusivities  of  the  two  components,  diffusive  plumes  grow 
across  the  interface  in  the  vertical  direction,  and 
transport  both  the  heat  and  salt  between  the  two  layers .  In 


3 


this  case,  the  faster  diffusing  constituent  (heat)  is 
destabilizing  and  provides  the  release  of  potential  energy 
to  drive  the  motion,  and  the  slower  diffusing  constituent 
(salt)  is  stabilizing.  The  regions  favorable  for  the 
diffusive  convection  include  the  polar  and  subpolar  regions 
of  the  Arctic  and  the  Southern  Ocean.  Here,  the  melting  of 
ice  and  small  rates  of  evaporation  result  in  the  colder  and 
fresher  upper  layer,  thereby  resulting  in  conditions 
favorable  to  diffusive  convection.  A  numerical  example  is 
shown  in  Figure  2 . 
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Figure  2.  Diffusive  convection  example.  Two  layer  system 

with  the  depiction  of  heat  (left)  and  salinity  (right) . 
The  high  values  of  T  and  S  are  shown  in  red;  low  values 

are  represented  in  blue. 


The  diffusive  convection,  in  turn,  takes  two  distinct 

forms  depending  on  the  background  stratification: 

oscillatory  convection  in  smooth  gradients  and  diffusive 

layering.  The  schematic  diagram  in  Figure  3  describes 

oscillatory  diffusive  convection.  In  this  model,  a  water 

parcel  is  slightly  displaced  upward  in  a  diffusively 

favorable  environment  (a  cold,  fresh  layer  overlays  a  warmer 

and  saltier  one)  .  As  the  particle  travels  upwards,  it  will 

rapidly  loose  heat  but  maintain  salinity.  However,  due  to 

the  decrease  of  the  background  salinity  with  depth,  the 
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parcel  becomes  heavier  than  its  surroundings,  and  the 
buoyancy  force  will  drive  it  back  towards  the  initial 
position.  The  parcel  becomes  heavier  from  the  initial  state 
not  only  because  of  the  change  in  the  density  contrast  with 
its  surroundings,  but  also  because  of  its  heat  loss.  As  a 
result,  the  restoring  buoyancy  force  which  drives  the  parcel 
to  its  initial  position  is  now  greater  than  the  initial 
upward  buoyancy  force.  The  final  result  is  a  gain  of  energy 
by  the  parcel  and  subsequent  oscillatory  motion  during  which 
the  parcel  overshoots  its  original  position  of  neutral 
equilibrium.  An  analogous  process  takes  place  as  the  parcel 
travels  downward  below  its  initial  position.  It  now  gains 
heat,  but  maintains  salinity,  and  since  the  parcel  is 
fresher  and  lighter,  it  will  experience  an  upward  buoyancy 
force  and  overshoot  its  initial  position.  The  whole  cycle 
will  be  repeated  time  and  again,  with  an  increase  of  energy 
within  the  parcel  and  the  growth  of  oscillations  over  time. 
This  instability  pattern  is  sometimes  referred  to  as  the 
"over-stable"  mode. 
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Cold,  Fresh 


Warm,  Salty 


Figure  3.  Schematic  diagram  of  oscillatory  diffusive 

convection . 

More  common  in  the  ocean  is  the  second  mode  of 
diffusive  convection,  known  as  the  diffusive  layering.  This 
mode  is  characterized  by  well  mixed  layers  separated  by 
sharp,  diffusive  interfaces.  Figure  4  represents  a  situation 
where  a  cold  and  fresh  layer  of  water  rests  on  the  top  of 
the  denser  warm  and  salty  layer.  In  this  configuration,  the 
heat  will  diffuse  faster  than  salt,  resulting  in  the 
downward  density  flux  across  the  interface.  The  region 
directly  below  (above)  the  interface  becomes  denser 
(lighter) .  As  a  result,  the  top-heavy  convection  is 
maintained  on  both  sides  of  the  interface. 
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Cold,  Fresh 

Thin  interface 

Ah 

Warm,  Salty 

Figure  4.  Schematic  diagram  of  diffusive  layering. 


B.  REGIONS  SUCEPTABLE  TO  DOUBLE  DIFFUSION 


In  order  to  determine  which  part  of  the  world' s  oceans 
will  be  favorable  for  double  diffusion,  two  different 
parameters  can  be  used.  Turner  (1965)  described  the 
dependence  of  the  relative  strength  of  double  diffusion  on 
the  density  ratio  as 


_  aTz 

PSZ 


(1) 


where  a  is  the  coefficient  of  thermal  expansion,  /?  is  the 
coefficient  of  saline  contraction,  and  Tzand  Sz  are  the 
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vertical  temperature  and  salinity  gradients,  respectively. 
The  definition  of  density  ratio  (1)  is  used  for  the  salt 
finger  regime  of  double  diffusion;  the  diffusive  convection 
is  described  by  the  inverse  of  (1) .  The  most  intense  double 
diffusion  for  both  regimes  is  realized  when  Rp  approaches 
1.  There  are  associated  difficulties  with  mapping  world 
oceans'  favorability  to  double  diffusion  based  on  Rp .  You 
(2002)  points  out  that  most  ocean  regions  often  show 
extremely  large  positive  or  negative  values  of  Rp ,  making 
the  interpretation  difficult. 


the 


To  simplify  interpretation,  Ruddick 
Turner  angle,  Tu,  named  in  honor  of  J. 


Tu (deg)  =  tan  1 


f  dT  ndS  dT  n 
a- — P—,a  —  +  p 


dz 


dz  dz 


(1983) 

Stewart 

as ^ 
dz  J 


presented 
Turner : 

(2) 


Rp  and  Tu  can  be  related  by 

Rp  =  -tan(lw  +  45)  .  (3) 

The  physical  meaning  of  the  Turner  angle  is  described 
in  Figure  5.  To  summarize,  when  -90  deg  <Tu<-45  deg  (or 
0<R  <1,  as  defined  in  equation  (1),  diffusive  convection  is 

possible;  when  -45  deg  <Tu<45  deg  (or  0<R/,<  ±  oo  )  ,  the  water 

column  is  stably  stratified  with  respect  to  T  and  S,  and 
double  diffusion  is  not  possible;  and  when  45  deg  <Tu<90 
(l<Rp<±oo)  the  salt  finger  regime  can  be  expected. 

You  (2002)  used  the  1994  Levitus  climatological  atlas 
to  calculate  Tu  for  the  world  oceans.  The  three  major  ocean 
basins  are  displayed  in  Figures  6-8.  He  notes  that  about  44% 
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of  the  oceans  are  favorable  for  double  diffusion,  of  which 
30%  takes  the  form  of  salt  fingering,  and  14%  is  diffusively 
favorable . 


Tu  *  45° 


Figure  5.  Illustration  of  Turner  angle  Tu  (degrees ) with 
corresponding  Revalue  indicated.  Adopted  from  You 

(2002)  . 
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Depth  (m) 


Figure  6.  3D  plots  of  (a)  Tu :  the  color  bar  marks  diffusive 
convection (-90  deg<Tu<-45  deg),  doubly  stable  (-45 
deg<Tu<45  deg),  and  salt-fingering (45  deg<Tu<90 
deg), (b) diffusive  convection  only  and  (c) salt-fingering 
only  in  the  western  (left  panel)  and  eastern  (right 
panel)  Altantic  Ocean.  Adopted  from  You  (2002)  . 
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Figure  7 


Same  as  Figure  4,  but  for  Indian  Ocean 
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Figure  8 


Same  as  Figure  5,  but  for  the  Pacific  Ocean. 
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C .  THERMOHALINE  STAIRCASES 

One  of  the  most  fascinating  aspects  of  double-diffusive 
convection  is  related  to  its  ability  to  form  stepped 
structures.  These  so  called  thermohaline  staircases  consist 
of  mixed  layers  which  are  separated  by  thin  stratified 
interfaces.  There  have  been  several  theories  proposed  to 
explain  the  origin  of  theromohaline  staircases.  One  of  the 
earliest  is  the  collective  instability  mechanism  (Stern, 
1969) ,  which  attributes  the  staircase  formation  to  the 
interaction  between  waves  and  double-diffusion,  resulting  in 
a  system  of  sheets  and  layers;  a  relatively  recent 
hypothesis  that  steps  in  temperature  and  salinity  profiles 
represent  intrusions  evolving  into  a  staircase  (Merryfield, 
2000) ;  and  the  latest  theory  proposes  that  layers  form  as  a 
result  of  the  instability  of  the  flux  gradient  laws  -  the 
instability  which  manifests  itself  in  the  form  of  growing, 
horizontally  uniform  perturbations  transforming  the  basic 
gradient  into  a  well-defined  thermohaline  staircase  (Radko, 
2003,  2005) . 

Salt  finger  staircases  are  normally  observed  when  Rp 

is  less  than  2;  however,  staircases  associated  with 
diffusive-convection  have  been  observed  over  a  wider  range 
of  l<Rp<10.  Recent  field  data  from  the  Atlantic  Ocean  have 

shown  that  vertical  fluxes  within  thermohaline  staircases 
are  significantly  (by  an  order  of  magnitude)  enhanced  over 
non-staircase  regions  with  similar  stratifications  (Schmitt 
et  al .  ,  2005)  .  Thus,  it  is  vital  to  improve  our 
understanding  of  thermohaline  staircase  formation,  evolution 
and  the  large  scale  consequences.  In  this  study,  it  is 
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argued  that  staircases  could  have  an  impact  on  the  oceanic 
heat  budget  in  the  Arctic  Ocean. 

D.  EQUILIBRIUM  STRUCTURE  OF  DIFFUSIVE  STAIRCASES 


Another  interesting  feature  of  diffusive  staircases  is 
related  to  their  evolution  over  time.  The  layers  which  form 
initially  are  unstable  and  thin.  They  merge  continuously, 
which  results  in  increasing  step  height  until  an  equilibrium 
thickness  is  reached.  This  merging  of  layers  in  diffusive 
staircases  is  attributed  to  the  secondary  instabilities  of 
steady-state  staircases  (Radko,  2005)  .  Two  merging 
mechanisms  are  possible:  either  layers  merge  when  interfaces 
drift  and  collide  (Balmforth  et  al .  ,  1998),  or  the  steps 
with  small  buoyancy  variations  become  even  weaker  and 
eventually  vanish  (Merryfield,  2000).  Radko  (2005)  observed 
that  formation  and  equilibration  of  the  staircase  increased 
the  diffusive  fluxes  by  almost  an  order  of  magnitude. 


One  of  the  aims  of  this  study  is  the  evaluation  of 
diffusive  fluxes;  therefore,  it  is  necessary  to  introduce 
the  flux  law  for  diffusive  convection.  The  so-called  4/3 
flux  law  for  diffusive  convection  was  formulated  by  Turner 
(1965)  .  The  4/3  flux  law  relates  temperature  flux  FT  (in 
degrees  C  m  s.-1)  and  salinity  flux  Fs  (in  m  s._1)  to  the 
temperature  and  salinity  steps  across  the  diffusive 
interface  (  AT  and  AS  )  as  follows: 


aFT  =  C(g/r3/v)1/3(aAT)4/3  (4) 

/?Fs=aF \Y  (5) 


where  C (R)  and  y{R  ),  kT  is  the  molecular  diffusivity  of 

P  P 

heat,  v  is  the  kinematic  viscosity,  g  is  the  gravitational 
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acceleration,  a  is  the  coefficient  of  thermal  expansion,  ft 
is  the  coefficient  of  haline  contraction,  y  is  the  flux 
ratio,  and  C  is  a  constant  determined  from  experimental  data 
or  oceanic  observations  assuming  that  the  4/3  flux  law  is 
valid.  In  (4)  and  (5),  dimensional  fluxes  are  used.  The  4/3 
flux  law  is  therefore  independent  of  stair  thickness,  but 
dependent  on  the  density  ratio.  Dimensional  considerations 
(Turner,  1973)  suggest  that  both  y  and  C  are  uniquely 

determined  by  the  density  ratio,  R  . 
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II .  NUMERICAL  SIMULATION 


A.  MODEL  DESCRIPTION 


In  order  to  study  oscillatory  diffusive  convection,  two 
dimensional  numerical  simulations  were  performed  on  the 
Naval  Postgraduate  School  High  Performance  Computing  Center 
Linux  clusters.  Following  Radko  &  Stern  (1999),  the 
temperature  and  salinity  fields  were  decomposed  into  the 

linear  basic  state  ( T,  S )  and  perturbations  ( S, T  )  of  it.  For 
computational  convenience,  Boussinesq  equations  were  non- 
dimensionalized.  Therefore,  d=(kTV'/ga71-)1/4  represents  the  unit 
of  length,  kT/d  is  the  velocity  scale,  d2/kT  represents  the 
time  scale,  kT  is  the  molecular  diffusivity  of  heat,  v  is 
the  kinematic  viscosity,  g  is  the  gravitational 
acceleration,  a  is  the  coefficient  of  thermal  expansion  and 

Tz is  the  temperature  gradient.  Wilson  (2007)  computed 

Tz  =0.015  based  on  the  average  1.5  deg  C  decrease  in  T  over 
100  m  depth  in  the  diffusive  layer  in  the  Beaufort  Gyre. 
Similarly,  aTzd  is  the  scale  for  both  temperature  and 
salinity  perturbations,  resulting  in 


—  —  =  -Vp  +  V2v  +  (T-  S)k, 
Pr  dt 

V  •  v  =  0, 


— -ir  =  v2r, 

dt 

dS  1 

- w  =  rv  S, 

dt  R „ 


(6) 
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is  the 


where  Pr=WkT  is  the  Prandtl  number,  r=ks/kT 
molecular  diffusivity  ratio  (Lewis  number) ,  R  is  the 

density  ratio  for  the  basic  uniform  S  -  T  gradient,  and  v, 
T,  S  are  non-dimensional  velocity,  temperature,  and 
salinity.  For  parameters  of  Arctic  staircases,  d~ 1  cm, 

t~ 686  seconds. 

The  conversion  formula  for  temperature  flux  is 

^Tdim  —  ^Trion-dim  XydT/ dz  (7) 

and  the  corresponding  heat  flux  is 

FHdim  (W  m'2)  =  FTdim  (m/s  K)CP  p  (8) 

where  Cp  is  the  specific  heat  for  sea  water  and  p  is  the 
density  of  sea  water.  Substituting  (7)  into  (8)  yields  the 
conversion  factor  of  0.0084. 

Unless  specified  otherwise,  the  diffusivity  ratio  of 
T  =  1/10  was  used,  which  is  an  order  of  magnitude  greater 
than  the  true  value  for  heat/salt  ( z  =  1/100 )  .  This  greater 
value  was  selected  purely  for  numerical  convenience. 
Numerical  calculations  with  a  true  diffusivity  ratio  would 
require  resolving  the  extremely  fine  salt  dissipation  scale, 
which  is  by  far  less  than  the  characteristic  size  of  the 
diffusive  layers,  and  the  calculations  would  be 
computationally  prohibitive.  As  was  pointed  out  in  Stern  et 
al . ,  (2001),  such  a  modification  is  not  expected  to  alter  the 
fundamental  physics  and  dynamics  of  the  diffusive 
convection.  Furthermore,  a  simple  extrapolation  procedure  is 
be  used  to  evaluate  the  heat/salt  (r  =  0.01)  fluxes  from  r  =  0.1 
numerical  simulations. 
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A  fully  dealiased  pseudospectral  method  described  and 
first  used  in  Stern  &  Radko(1998)  is  employed  to  obtain  a 
two  dimensional  numerical  solution  for  the  system  of 
equations  (6)  above.  The  solution  is  doubly-periodic  in  x 
and  z.  All  numerical  algorithms  were  coded  in  Fortran. 

B.  PRELIMINARY  CALCULATIONS 

In  order  to  gain  confidence  in  the  numerical 
calculations,  and  to  better  understand  the  dynamics  of  the 
diffusive  convection,  a  simple  experiment  was  carried  out 
initially.  The  experiment  was  initiated  from  a  random 
computer-generated  distribution  of  T  and  S  and  is  shown  in 
Figures  9  and  10.  This  experiment  was  performed  using  a 
uniform  grid  of  512x1024  nodes  and  non-dimensional  grid 
spacing  delta  x=y=0.75  which  corresponds  to  a  vertical  scale 
of  7.5  m,  and  =  1.1.  The  first  stage  of  the  experiment 

(t<50)  is  characterized  by  the  development  of  small  scale 
diffusive  plumes  presented  in  Figure  9.  At  t=500,  the  field 
of  diffusive  plumes  arranged  itself  into  a  system  of  well 
defined  horizontal  steps.  Next,  a  series  of  merging  events 
occurred:  between  t=  500  and  t=600  the  number  of  steps 
decreased  from  9  to  6  as  shown  in  Figures  9  and  10.  Six 
steps  merged  into  three  in  less  than  100  time  units,  or 
approximately  19  hours.  The  experiment  reached  the 
equilibrium  at  approximately  t=800  when  the  staircase 
evolved  into  two  final  steps,  as  shown  in  Figure  11.  The 
equilibrium  step  height  in  the  final  stage  was  3.75  m.  The 
experiment  was  continued  until  t=950;  however,  no  merging 
took  place  during  this  period. 
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The  time  series  of  fluxes  for  this  experiment,  which 
are  computed  over  the  whole  computational  domain,  is  shown 
in  Figure  12.  It  confirms  that  the  diffusive  fluxes  increase 
with  each  merging  event,  and  once  the  equilibrium  thickness 
of  the  steps  is  reached,  the  fluxes  decrease  by  about  half 
of  their  value  and,  subsequently,  remain  roughly  constant. 

This  preliminary  numerical  experiment  indicates: 

(i)  sensitivity  of  fluxes  to  the  step  size 

(ii)  the  tendency  of  the  weak  steps  to  merge 

The  mechanics  of  merging  in  a  staircase  and  associated 
fluxes  are  now  analyzed  in  greater  detail. 
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Figure  9.  Formation  of  layers  in  a  double-diffusive 

numerical  experiment.  The  left  figure  shows  a  well 
developed  field  of  two  dimensional  diffusive  plumes 
after  t=50.  Note  the  well  defined  horizontal  layers  in 
the  right  figure  after  t=500.  The  red  color  corresponds 
to  high  values  of  T;  the  low  values  are  shown  in  blue. 


21 


Figure  10.  Continuation  of  the  numerical  experiment  from 

Figure  9  demonstrates  merging  of  diffusive  layers.  Note 
six  layers  in  the  left  figure  at  t=600  merging  into 
four  layers  at  t=700. 
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Figure  11.  Continuation  of  the  numerical  experiment  from 

Figures  9  and  10.  The  equilibrium  state  shows  2  layers. 


23 


Figure  12.  Time  series  of  diffusive  fluxes  in  the  numerical 
experiment  initiated  from  random  noise.  The  blue  curve 
indicates  heat  flux,  the  green  one  salinity  flux.  The 
arrows  indicate  a  merging  event,  the  red  arrow  an 

equilibrium  state. 


C  ONE -STEP  FLUXES 

Turner  (1973)  extensively  discussed  the  diffusive 
instability.  One  of  the  most  important  properties  of  well 
formed  diffusive  interfaces  is  the  temperature  and  salinity 
flux  across  them.  Padman  (1994)  and  Muench  et  al  (1990) 
indicated  that  diffusive  convection  could  be  a  significant 
process  influencing  the  vertical  fluxes  of  heat  and  salt 
from  the  warm  and  salty  intermediate  water  upwards  to  the 
polar  mixed  layers  and  sea  ice.  The  basin-averaged  estimate 
of  heat  loss  from  the  Atlantic  layer  is  7  W  rrf2  (Aagaard  and 
Greisman,  1975).  Previous  studies  (Padman  et  al .  ,  1987) 


estimated  heat  flux  using  the  extrapolation  of  the 
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laboratory  flux  laws  to  the  ocean  and  reported  values  of 
0.02  to  0.1  W  m-2 .  However,  based  on  the  analysis  of  the 
Arctic  thermohaline  staircase  from  the  Beaufort  Gyre 
experiment,  Wilson  (2007)  reported  that  the  diffusive  fluxes 
could  be  an  order  of  magnitude  greater  than  the  diffusive 
fluxes  determined  from  previous  laboratory  experiments.  Use 
of  the  direct  numerical  simulations  opens  a  possibility  for 
resolving  the  controversy  with  regard  to  the  magnitude  of 
diffusive  fluxes. 

A  series  of  one-step  numerical  experiments  was 
performed  in  order  to  determine  the  magnitude  of  diffusive 
fluxes  numerically  and  compare  them  with  those  from  the 
laboratory  and  field  observations.  The  range  of  Rp  was  from 

1.1  to  5.0.  The  upper  range  of  the  values  of  R  is  more 
significant  for  the  comparison  with  the  Wilson  (2007)  Arctic 
Ocean  results,  since  the  values  of  R  in  the  thermohaline 

staircases  in  that  geographic  region  are  5  or  higher.  The 
vertical  scale  of  steps  was  varied  as  well,  in  order  to 
examine  dependencies  of  simulated  fluxes  and  structures  on 
H.  The  step  thicknesses  were  varied  from  0.75  to  3  m  and  the 
results  are  summarized  in  Table  1. 

As  expected,  the  fluxes  generally  decrease  with 
increasing  R  .  The  fluxes  also  generally  increase  with 
increasing  layer  thickness,  H,  except  for  a  few  values  in 
the  range  of  R  <1.9  and  H>3.0m.  However,  these  differences 

are  small  and  can  be  attributed  to  the  incomplete 
equilibration  in  some  runs. 
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Layer 

Thickness 

(m) 

R, 

1.1 

1.3 

1.5 

1.7 

1.9 

2.1 

2.3 

3.0 

4.0 

5.0 

0.75 

62 . 8 

28 . 9 

20.9 

13.1 

9.13 

6.90 

4 . 94 

3.08 

1 . 95 

1.24 

1.5 

155 

58.2 

43.1 

33.2 

24 . 9 

20.5 

20.8 

11 . 8 

9.27 

7.76 

3.0 

45.9 

49.2 

31.3 

30.3 

26.5 

18.0 

19.0 

14 . 1 

11.3 

10.3 

Table  1.  Diffusive  heat  flux  (non-dimensional)  as  the 

function  of  density  ratio  and  layer  thickness. 

D  MERGING  EVENTS 

Next,  in  order  to  better  understand  the  mechanics  of 
merging  in  a  staircase  (see  Figures  9,  10) ,  several 

numerical  experiments  were  performed  to  isolate  the  merging 
process.  The  initial  setup  consisted  of  two  almost  identical 
steps;  the  focus  was  on  the  evolutionary  pattern  and  the 
time  scale  of  expected  merging  events.  All  experiments  were 
initiated  in  the  mesh  grid  of  512  x  512  points  and  delta 
x=y=0 . 3  which  corresponded  to  a  vertical  scale  of  1.5  m.  The 


value  of 

R, 

was  varied  from 

1.2  to 

1.5.  For 

each 

experiment. 

the 

value  of  was 

slightly 

higher  for  one 

interface 

and 

slightly  lower 

for  the  other. 

The 

representative  numerical  experiment 

in  this 

series  is 

shown 

in  Figures  13-15.  After  approximately  t=100,  the  erosion  of 
one  of  the  steps  is  noticeable  in  Figure  13.  The  process 
continues  until  the  step  is  completely  eroded  at  t=400. 

Figures  14  and  15  illustrate  a  gradual  weakening  of  the 
interface  until  its  complete  disappearance  at  t=400. 
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indicating  that  a  merging  event  took  place  between  t=300  and 
t=400.  Visual  inspection  of  the  evolutionary  pattern 
indicates  that  the  interfaces  remains  stationary  and  the 
weaker  one  is  gradually  eroded  until  it  completely 
disappears . 

Figure  16  represents  a  time  series  of  numerical  fluxes 
for  the  same  two-step  experiment.  Note  the  increase  in  the 
fluxes  indicating  the  merging  of  two  steps  into  one.  The 
jump  in  fluxes  occurs  at  t=350,  indicating  the  exact  time 
the  merging  event  took  place.  This  coincides  with  findings 
from  Figures  13-15. 

The  experiments  are  summarized  in  Table  2.  The  results 
confirm  that  with  increasing  R  ,  the  diffusive  fluxes 

decrease.  Each  experiment  exhibits  a  rough  doubling  of 
fluxes  once  merging  takes  place.  The  time  scale  of  merging 
increases  with  increasing  R  .  The  numerical  results  will  be 

explored  further  once  the  theoretical  framework  for 
understanding  merging  events  -  the  merging  theorem  -  is 
developed  in  Chapter  III. 
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Figure  13.  Numerical  experiment  showing  merging  evolution  of 
two  steps  into  one.  The  plots  show  the  total  averaged 
density  for  the  simulation,  from  left  to  right  at 
t=25, 100, 200, 300, 350, 400 .  The  initial  step  height  in 
the  experiment  is  0.75  meters. 
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Figure  14.  Evolution  of  merging  event  in  a  numerical 

experiment  of  a  two-step  system.  The  top  figure  shows  a 
sharp  interface  and  clearly  formed  diffusive  plumes. 
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Figure  15.  Continuation  of  merging  experiment.  The  top  figure 
shows  weak  interface  at  t=300.  In  the  bottom  figure, 
the  middle  interface  is  no  longer  present  at  t=400. 
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Figure  16.  Time  series  of  diffusive  fluxes  in  a  two-step 
numerical  experiment.  The  blue  curve  represents  heat 
flux;  the  green  one  represents  salinity. 


Ft  before 

Ft  after 

Time 

merging 

merging 

(non-dim. ) 

1.2 

48 

119 

130 

1.3 

41 

83 

350 

1.4 

31 

70 

950 

Table  2.  Summary  of  two-step  numerical  experiments. 
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Ill .  MERGING  THEOREM 


A.  INTRODUCTION 

The  existence  of  thermohaline  staircases  in  the  ocean 
has  been  well  established  (Kelley,  1984;  Schmitt,  1994). 
They  frequently  occur  in  the  ocean  regions  where  the 
conditions  for  double  diffusive  instability  are  favorable. 
The  mechanics  proposed  for  their  origin  and  evolution  have 
been  discussed  in  the  Introduction.  Here,  the  mechanism  for 
the  merging  of  such  staircases  will  be  explored  and 
subsequently  tested  in  a  series  of  numerical  experiments. 

Two  mechanisms  for  the  merging  of  staircases  have  been 
studied:  either  layers  merge  when  interfaces  drift  and 
collide,  or  the  interfaces  will  erode  over  time  without 
exhibiting  any  vertical  motion.  The  first  mechanism  has  been 
proposed  by  Balmforth  et  al  (1998)  who  formulated  a  one¬ 
dimensional  model  of  the  mechanically  forced  turbulence  in  a 
stratified  fluid.  Following  Radko  (2007),  this  mechanism 
will  be  referred  to  as  H-merger.  The  second  mechanism,  which 
was  described  by  Merryfield  (2000) ,  who  modeled  the  up 
gradient  buoyancy  flux  in  a  double-diffusion  favorable 
fluid.  According  to  his  model,  the  steps  with  small  buoyancy 
variations  weaken  even  further  and  eventually  vanish.  Again, 
following  Radko  (2007),  this  mechanism  will  be  referred  to 
as  B-merger. 

This  chapter  is  organized  as  follows:  First,  the 
merging  theorem  for  the  two-component  H-  and  B-  merger  will 
be  formulated.  Next,  a  physical  explanation  of  B-  merger 
will  be  presented,  and  the  relationship  between  the  density 
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will  be 


ratio,  R  ,  and  the  diffusivity  ratio,  y, 

demonstrated.  Lastly,  the  relationship  between  the  mechanics 
of  merging  and  the  variation  in  the  diffusivity  ratio  will 
be  explored. 

The  starting  point  for  the  derivation  of  the  two 

component  merging  theorem  is  the  two  component  buoyancy 

conservation  equation  where  T(z,t)  and  S(z,t)  are  heat  and 

salinity,  respectively.  FT  and  Fs  are  heat  and  salinity 

fluxes,  respectively. 

dT  _  dFs ' 

dt  dz  /  g  \ 

dS  dFT  I 
dt  dz  , 
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Figure  17.  Schematic  representation  of  the  stability  analysis 
for  a  periodic  system  of  layers  and  interfaces,  (a) 
Basic  state  with  equal  steps,  (b)  Perturbed  state  with 
the  slightly  increased  buoyancy  jumps  at  even 
interfaces  and  with  the  correspondingly  decreased  jumps 

at  odd  interfaces. 


Following  Radko' s  (2007)  one-component  theory,  we 
express  the  heat  (T)  and  salinity  (S)  variations  across  a 
two  step  system  in  Figure  17(b)  as  follows: 

T12~T01  =T-S 
Tn-Tn=T  +  8 

\  (10) 

S\2  ~Sm  =S-e 

S23-Sl2=S  +  s,  (8,s) « (T,S) 
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where  T  =  (dT  / dz)H  and  S  =  (dS  / dz)H  are  the  temperature  and 


salinity  variation  across  one  step  in  Figure  17  (a) , 
respectively.  Symbols  5  and  s  represent  a  small 
perturbation  to  T  and  S,  respectively.  Note  that  the 
thicknesses  of  two  layers  in  Figure  17  (b)  are  not 

necessarily  equal: 

Hl2=H-h 

H23  =  H  +  h,  h  «  H 


where  Hnn+i  represents  the  distance  between  the  two 
interfaces  z  =  and  z  =  zn+l.  For  convenience,  the  buoyancy 
equation  is  expressed  in  the  integral  form,  which  is 
accomplished  by  integrating  (9)  over  the  interval  [ zbot,  ztop]  : 


7  dT  z7  dFT 

J  J-dZ  =  J  —dZ  =  F  (Z">p)-F  (hot) 


zbot 


dt 


’■lop  :,op 

l-y'z=lv^z=FS(zb-F7w,> 


dt  :  dt 

zbot  zbot 


(12: 


where  two  levels  ztop(t)  and  zbot(t)  can  vary  in  time 
Utilizing  the  identity 


dy  y  QT  d  d 

J  Tdz=\  — ^+—n%) — 


dt 


zbot 


zbot 


dt 


dt 


J  J  Qf  V  ,OP>  rtf  V  bot> 


dt 


zbot 


zbot 


dt 


dt 


(i3: 


the  equation  (13)  is  re-written  to  a  new  form 


ij Tiz  - 

Zbot 

7  Ztop  J_  r 

zbot 


(14) 
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Next,  the  integral  relation  in  (14)  is  applied  to  the 


individual  layers  [zn,zn+1  resulting  in 


d_ 

dt 


zn+ 1 

J  Tdz 


zn 


d_ 

dt 


zn+ 1 

J  Sdz 


Zn 


n+\  rri  _  rr  ■  tjT  _  tjT 

j ,  n+ 1  ” n+ 1 

at  at 


dZn+\  c  ^  n  C  i  /7s  _  J7S 

,  ^n+l  ,  'h  ^n+1 

dt  dt 


(15) 


where  Fi=F(zi)  .  The  greatest  contribution  to  the  integrals 
of  temperature  and  salinity  in  equation  (15)  comes  from  the 
interior  of  layers  where  the  gradients  of  temperature  and 
salinity  are  nearly  uniform,  and  therefore 


zn+ 1 

J  Tdz*T.^H. 

Zn 
Zn+ 1 

J  Sdz*S,„,H, 


(16) 


Applying  equations  (15)  and  (16) to  two  successive  layers 
[zn,zn+l]  where  n  =  1,2  yields 

ji(HnTn)  =  F2T-F,r 

(17) 

jt(Hi  A,)=d-d- 

Due  to  a  substantial  variation  in  temperature  and  salinity 
across  the  interfaces,  it  is  essential  to  define  the 
temperature  and  salinity  at  the  interface  Tn,Sn  .  The  center 
of  the  interface  is  defined  as  a  point  where  the  temperature 
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and  salinity  are  average  values  at  the  center  of  the  two 
neighboring  layers 


(18) 


A  similar  set  of  integral  equations  results  from  the 
temperature  and  salinity  budgets  for  areas  enclosed  by  the 
centers  of  the  neighboring  layers  [znn+1]  : 


where  Fm+l  =  F(znn+1);  (AT,  AS)L12  and  (AT,  AS)L23  are  the  temperature 
and  salinity  variations  across  the  interior  of  the  layers  in 
Figure  17  (b)  .  To  formulate  the  equations  for  evolution  of 
S,s,  and  h  in  time,  the  equations  (17)  and  (19)  are  utilized 
and  the  result  is  simplified  using  periodicity  conditions, 
yielding : 


H  —  +  2h ^ 
dt  dt 


^  =  -2(F2T-Flr) 


H  —  +  2h 
dt  dt 


-r  =  -2  (f2s-f?) 


Since  it  is  assumed  that  the  perturbations  are  relatively 
small  when  compared  to  the  basic  state,  (20)  is  linearized: 
We  neglect  the  nonlinear  terms  h(dT2/dt)  and  h(dS2/dt)  in 
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(20),  further  reducing (20)  to 


II 

'c 

'4-4) 

dt 

H[ 

K  2  1  / 

d 

—  £  = 

4-4) 

dt 

HX 

,  2  1  / 

the 

two 

equati 

(2i ; 


simplifying  using  the  periodicity  conditions,  the  following 
is  obtained: 


dh 

8Tlu  +  A  Tl22 

dh 

H  d 

dt 

4 

dt 

4  dt 

dh 

^Z12  +  8Sl  23 

dh 

H  d 

dt 

4 

dt 

4  dt 

-■ T ■ T 0 ^23  - ■ ATLn )  =  2(4  - 4 )  +  2 


JL  23 


JL12> 


Vs  . 
23 


12  > 


dZy 

dt 
dz , 


12 


dt 


8 , 


(22: 


Equation  (22)  is  linearized  about  the  steady  state  in  Figure 

dz 

17.  Since  dz„ldt  =  0  in  basic  steady  state,  terms  — —  f  and 
12  dt 


dz. 


dt 


■8  approach  zero  in  the  limit  of  a  weak  perturbation.  It 


is  reasonable  then  to  neglect  the  non-linear  term  2 (dzv  / dt)  . 
Realizing  that  (AT,AS)Ll2~(AT,AS)L23~(AT,AS)L,  the  h-equation  is 
reduced  to 


f 

T 

v 


( 

S 

v 


A7p 

2  J 

2  J 


dh  H  d 
dt  4  dt 


{8Tl  23 


dh 

dt 


tL£iAS 

4  dty  J 


L  23 


-Arm)  =  2  (F2p^) 
-A5m)  =  2  (Fi~FS). 


(23) 


Replacing  the  instantaneous  fluxes  in  layers  and 
interfaces  andF^  in  (21)  by  the  steady  one-step  fluxes 

Ft  and  Fs  (Radko  2007)  based  on  the  parameters  of  a  region 
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extending  between  the  centers  of  two  adjacent  layers 
(zn-in  andz„n-i) '  we  arrive  at: 


F1  »F  (Tn  —  T0l,Sn  —  SQX,zn  —  z01) 
F2  &  F  (T23  —Tt2,S23  -  S]2,z23  —  zi2) 
ij  »F  (Tn  —  Tai,Sn  —  Sm,zn  —  z0i) 
^2  ~  F  (^12  _  ^01’Z12  —  Z0l)’ 


(24) 


Next,  the  difference  between  the  temperature  and  salinity 
fluxes  at  the  adjacent  interfaces  is  reduced  to 


8Ft 

dPT ' 

F2t-Ft  =25 

+  2  £ 

Z  I 

df 

dS 

dFs 

dFs 

F2  -F,s  =25 

+  2  £ 

Z  1 

df 

dS  , 

and  (25)  becomes 

d  45  dFT  4s  dFT ' 
dt  ~  H  df  H  8S 
d_  45  dFs  4s  dFs 
dt£~  H  df  H  DS  , 


1 .  H  -  Merger 

The  fluxes  at  the  center  of  each  layer  Fnn+l  are 
approximated  by  their  corresponding  one  step  layer  based  on 
the  thickness  and  buoyancy  variation  across  this  layer 

FZ=Ft(T2-T1,S2-Suz2-z1)' 

F^=Ft{T3-T2,S3-S2,z3-z2)\  (27) 

Ff2  =  Fs(T2-Tl,S2-S1,z2-zl) 
F2s3=Fs(T3-T2,S3-S2,z3-z2 l 
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J 


( 

< 

*- 

V 


dh 

H  d 

dt 

4  dt 

dh 

H  d 

dt 

4  dt 

(A  Tl23~ 


2 

\S 
2  ) 


AT,„-AT,„  h 


ATL12)  =  4h i- 

il2;  a// 


a# 


(28: 


lL  23 


Z12 


a// 

aA.s, 


ASL23-ASLU=^F2h 
oH 


(29) 


and  finally,  substitution  of  the  normal  modes 
(S,£,h)  =  (S0,£0,h0)exp(At)  in  equation  (28)  results  in  the 
eiegnevalue  equation  for  H-growth  rates: 


Ah  — 


dFT 

dH 


8FS 

8H 


A Tl  H  8 AT,  AS,  H  8 AS,  ' 

1 - ~ —  o  — 


oo : 


2  2  8H  2  2  8H 

This  equation  indicates  that  if  fluxes  are  insensitive 

to  the  layer  thickness,  then  AH  =  0  and  the  H-merger  events 
cannot  occur.  This  observation  may  rationalize  the  absence 
of  H-mergers  in  the  numerical  simulations  (Figures  9-11) . 

2 .  B  -  Merger 

Substitution  of  the  normal  modes  (8,£,h)  =  (S0,£0,h0  )exp(/?f) 
in  (26)  results  in  the  eigenvalue  equation  for  growth  rates: 


ab+— 

B  H 


8F  8F 
+  - 


8T  8S 


5  16 

A„  H - T- 

*  H2 


8Ft  8Fs  8Fs  8F 7 


8T  8S  8T  8S 


=  0 


(3i : 


Since  the  free  coefficient  is  expected  to  be  small,  the 
positive  root  can  be  approximated  by 
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(32) 


dFs  dFT  dFT  dFs 

_  4  df  ds  df  dS 
H  dFT  ^  dFs 

1)F+1)F 


3 .  Specific  Solutions 


Suppose  that  Turners'  4/3  flux  laws  are  relevant  for 
diffusive  convection  and  can  be  adopted  for  the  analysis  of 
numerical  solutions.  (The  numerical  results  generally 
support  this  view.)  The  4/3  flux  law  from  Turner  (1973)  has 
a  form 

4 

FT=C(Rp)AT>  (33) 

Fs=r(Rp)FT 

Substituting  (33)  into  (32)  yields 

C2  5/ 

x  jmrr - ^ -  .  ,34, 

H  3C'Rp-4C-3C'r-3Cy' 


The  last  two  terms  in  the  denominator  of  (34)  are  very 
small,  and  therefore  can  be  neglected.  This  simplifies  (34) 
to 


C2  dy 

A  _  16(A T)m  dRp 
B  H  3C'Rp-4C 


(35) 


Since  C\R  )<  0,  AB  is  positive 
decreases  with  .  This 

staircases  are  susceptible  to 


only  if  the  flux  ratio  y(Rp) 
implies  that  thermohaline 

dy 


B-merger  as  long  as 


6R 


<0 
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All  of  our  numerical  simulations  summarized  in  figure  18 
indicate  that  this  condition  is  satisfied  and  thereby  (35) 
rationalizes  numerous  B-merger  events  as  discussed  in 
Chapter  II. 

On  the  other  hand,  since  fluxes  in  (33)  are  independent 
of  step  height,  AH  in  (30)  is  identically  zero,  which 

rationalizes  the  absence  of  H-merging  events  in  our 
preliminary  numerical  simulations. 


Figure  18.  Dependence  of  A  on  R  from  one  step  numerical 

experiments . 


43 


PHYSICAL  EXPLANATION 


B  . 


S(z) 

CO 

II 

Pi 

a¥T 

p  a  AT 

R 


p2 


Fs2<FT2 

FSi«FTi 


Interface  2 


Interface  1 


Figure  19.  Illustration  of  B-merger  and  y  dependence  on  R  . 

Depicted  is  a  two-step  staircase  where  FT  and  Fs  are 
heat  and  salinity  fluxes,  respectively.  The  dashed  line 
represents  erosion  of  the  step  due  to  the  greater  flux 
of  salt  across  interface  1  vice  2. 


Figure  19  attempts  to  explain  the  physical  processes 
involved  in  B-merger  of  the  thermohaline  staircase. 

The  starting  point  is  a  two-step  thermohaline 
staircase.  In  Figure  18,  two  steps  are  almost  identical,  and 
the  density  ratio  at  interface  1  ( Rpl  )  is  only  slightly 

greater  than  the  density  ratio  at  interface  2(R2). 
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Diffusive  convection  is  soon  established  due  to  the  fresh 
and  cold  layer  overlaying  the  warm  and  salty  layer, 
resulting  in  heat  and  salinity  fluxes,  FT  and  Fs,  across  the 
interfaces.  Now  it  is  necessary  to  introduce  the  idea  of 
dependence  of  the  density  on  the  flux  ratio,  y  .  As  Figure 
18  illustrates,  the  flux  ratio  decreases  with  the  increasing 
density  ratio.  Since  Rpl>  Rpl  ,  it  follows  that  the  interface 
1  has  a  slightly  lower  value  of  y  than  the  top  interface. 
Consequently,  the  salinity  flux  across  the  interface  1  is 
smaller  when  compared  to  the  interface  2.  The  net  result  is 
the  loss  of  salinity  from  the  bottom  step,  and  subsequent 
erosion  of  the  bottom  step.  As  salinity  flux  continues 
across  the  interfaces,  it  decreases  across  the  interface  1 
and  increases  across  the  interface  2,  thus  causing  an 
increase  of  AS  across  the  interface  1  and  decrease  of  AS 
across  the  interface  2.  This  in  turn  causes  Rp  across  the 

interface  1  to  become  even  greater,  and  across  the 

interface  2  becomes  even  smaller.  As  salinity  flux  continues 
to  increase  across  the  interface  2,  so  does  the  value  of  y 
for  the  same  interface.  Similarly,  the  salinity  flux  for  the 
interface  1  continues  to  decrease,  and  so  does  the  value  of 
y  for  that  interface. 

This  positive  feedback  leads  to  the  complete  erosion  of 
the  weaker  interface  2.  The  proposed  mechanism  is  consistent 
with  both  merging  theorem  (35)  and  with  the  earlier  (Chapter 
II)  numerical  simulations. 
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IV.  THEORY-BASED  ANALYSIS  OF  THE  DIRECT  NUMERICAL 

SIMULATIONS 


A.  TURNER'S  FLUX  LAW  MODEL 

1.  Verification  of  the  4/3  Flux  Law 

The  4/3  flux  law,  first  introduced  in  Chapter  I,  is 
used  in  calculating  of  the  oceanic  diffusive  fluxes  since 
their  direct  measurements  are  rarely  available.  Therefore, 
estimates  are  made  from  measuring  AT  and  calculating  C 
empirically,  then  substituting  into  (4) .  If  the  4/3  flux  law 
applies  to  our  numerical  calculations,  then  the  quantity 
F 

— rdim4  showed  in  Table  3  should  depend  largely  on  Rp  but  not 
«„„)> 
on  H . 

To  test  the  4/3  flux  law,  the  results  from  the  series 
of  one  step  numerical  experiments  introduced  earlier  were 
analyzed.  The  fully  equilibrated  results  of  numerical 
simulations  were  used,  and  the  averaged  heat  flux  was 
computed  based  only  on  the  second  half  of  the  experiment. 
The  goal  was  to  exclude  the  large  initial  fluctuations  of 
diffusive  fluxes  and  use  only  the  equilibrated  results.  The 
typical  non-dimensional  fluxes  from  a  one-step  numerical 
experiment  are  displayed  in  Figure  20.  Note  the  initial 
spike  in  the  flux  values. 

The  magnitudes  of  non-dimensional  fluxes  obtained  from 
each  experiment  were  first  dimensionalized  using  (7)  and 
(8)  .  Since  there  is  a  linear  relationship  between  AT  and 
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dT 

A H ,  both  quantities  are  equivalent  ( -  is  constant  in  our 

dz 

model)  and  AH  (vertical  step  scale)  was  used  in  equation 
(4) .  The  results  are  summarized  in  Table  3. 


GO 

li_ 

X 


Figure  20.  Time  plot  of  non-dimensionalized  diffusive  fluxes 
from  a  numerical  experiment.  The  blue  curve  and  green 
curve  represent  heat  and  salinity  flux,  respectively. 
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H 

(m) 

1 . 1 

1.3 

1.5 

1.7 

1 . 9 

2 . 1 

2.3 

3.0 

4.0 

5.0 

0.75 

.774 

.357 

.258 

.161 

.113 

.085 

.061 

.038 

.024 

.015 

1.5 

.758 

.285 

.211 

.162 

.  122 

.101 

.103 

.057 

0.45 

.037 

3.0 

.544 

- 

.  157 

- 

- 

- 

.083 

.035 

- 

.023 

Table  3.  — rdim4  as  a  function  of  the  density  ratio,  Rp  and 

n,m)> 

the  layer  thickness  (m) .  Dashes  indicate  no  values  were 
obtained  due  to  the  computational  constrains. 

An  analysis  of  the  results  from  Table  3  reveals  that 
the  4/3  flux  law  works  reasonably  well  for  our  numerical 
experiments  as  the  values  of  heat  flux  depend  more  on  Rp 

and  less  on  H. 

2 .  Comparison  of  Numerical  Fluxes  with  Fluxes 
Computed  from  Oceanic  Data 

There  is  some  disagreement  between  the  magnitude  of 
diffusive  fluxes  that  are  computed  using  the  laboratory 
derived  flux  laws  and  the  fluxes  computed  using  actual 
oceanic  observation.  If  the  diffusive  heat  fluxes  for  the 
Arctic  are  on  the  order  of  magnitude  2  to  6  W  rrf2,  as 

calculated  by  Wilson  (2007),  they  could  play  a  significant 
role  in  the  Arctic  heat  budget  and  contribute  to  the  melting 
of  the  Arctic  ice  cap. 

Padman  and  Dillon  (1987)  reported  that  based  on  the 

data  collected  from  the  Beaufort  Sea  during  the  Arctic 

Internal  Wave  Experiment  (AIWEX)  in  March-April  1985,  the 
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heat  fluxes  estimated  from  laboratory-derived  flux  laws 
ranged  from  0.02  to  0.1  W  rrf2  which  is  one  to  two  order  of 
magnitude  smaller  than  the  alternative  heat  flux  estimates. 
Padman  (1994)  estimates  vertical  heat  fluxes  through 
thermohaline  steps  to  be  between  2  and  5  W  rrf2  using  flux 
laws  proposed  by  Kelley  (1990)  and  Marmorino  and 
Caldwell ( 1 97 6)  .  His  findings  are  derived  from  the  data 
collected  near  the  Yermak  Plateau  as  a  part  of  the 
Coordinated  Eastern  Arctic  Experiment (CEAREX)  in  1989.  Based 
on  the  analysis  of  the  Ice-Tethered  Profiler  data  obtained 
from  Woods  Hole  Oceanographic  Institute  and  the  Beaufort 
Gyre  Exploration  Program,  heat  fluxes  calculated  using 
Turner's  4/3  flux  laws  range  from  2  to  6  W  rrf2  (Wilson, 
2007)  . 

However,  as  Padman  and  Dillon  (1988)  point  out, 
laboratory  experiments  differ  in  a  number  of  ways  from  the 
actual  ocean.  For  example,  much  greater  fluxes  are  required 
to  minimize  wall  heat  losses,  and  there  is  a  possibility  of 
secondary  circulations  driven  by  the  tank  geometry  on  the 
laboratory  side;  additionally  there  is  the  potential  for 
significant  interaction  between  double  diffusion  and 
instabilities  driven  by  ambient  shear.  Also,  laboratory 
derived  flux  laws  (Turner,  1973)  may  not  be  the  ideal  tool 
to  use  in  estimating  diffusive  fluxes  found  in  oceanic 
staircases.  Therefore,  it  is  reasonable  to  question  the 
assumption  that  the  values  of  diffusive  heat  fluxes  obtained 
from  laboratory  experiments  can  be  extrapolated  to  the 
actual  oceanic  conditions. 

Another  interesting  angle  to  consider  is  the  value  of 
the  molecular  diffusivity  ratio.  As  mentioned  in  Chapter  II, 
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the  numerical  experiments  relevant  to  this  thesis  used  the 
molecular  diffusivity  ratio  r  =  1/10  .  As  it  was  explained 
earlier,  this  was  done  for  numerical  convenience  and  should 
not  significantly  affect  the  physics  and  dynamics  of  the 
diffusive  convection.  It  is,  however,  highly  desirable  to 
take  into  account  dependence  of  fluxes  on  r  ,  which  was  done 
as  follows. 

Takao  and  Narusawa  (1980)  conducted  a  study  with 
substances  of  different  molecular  diffusivity  ratio,  z  ,  and 
concluded  that  the  substance  with  the  lowest  value  of  r 
registered  the  highest  flux  ratio.  They  suggested  that  the 
heat  fluxes  across  the  diffusive  interface  should  scale  as 

r~0'71  .  Thus  the  value  of  the  diffusive  fluxes  would  have  to 
increase  by  a  factor  of  5.13  (scaling  the  value  of  r  =  0.r"'71 
in  the  numerical  experiments)  .  In  addition,  in  order  to 
obtain  equivalent  values,  each  numerical  experiment  was 
extrapolated,  using  the  4/3  flux  law,  to  the  step  height  of 
5  m,  the  average  step  height  of  the  Beaufort  Gyre  Arctic 
staircase  (Wilson,  2007)  Table  4  summarizes  the  results 
after  applying  this  correction. 
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H 

(m) 

1 . 1 

1.3 

1.5 

1.7 

1 . 9 

2 . 1 

2.3 

3.0 

4 . 0 

5.0 

0.75 

33.8 

15.6 

2.20 

11.3 

5.1 

3.7 

2.7 

2.7 

1.0 

0.9 

1.5 

33.3 

12.3 

1 . 80 

9.23 

5.3 

4.7 

4.5 

2.5 

1 . 6 

1 . 6 

3.0 

18 . 1 

- 

5.28 

- 

- 

- 

3.4 

1 . 9 

- 

0.9 

Table  4.  Heat  flux  (W  m  2)  ,  corrected  by  the  factor  of 

5.13,  as  the  function  of  the  density  ratio,  Rp  ,  and  the 

layer  thickness,  H  (m) .  Dashes  indicate  no  values  were 
obtained  due  to  the  computational  constraints. 

Comparison  of  the  results  from  Table  4  to  the  diffusive 
fluxes  determined  from  oceanic  observations  by  Wilson 
(2007),  reveals  that  both  numerical  and  oceanic  heat  fluxes 
are  in  the  same  order  of  magnitude.  This  raises  the  question 
of  using  the  higher  value  of  r  in  the  numerical  simulations 
of  diffusive  heat  fluxes.  It  is  reasonable  to  assume  that 
the  higher  value  of  r  adversely  affects  the  numerically 
determined  diffusive  heat  fluxes  and  warrants  future 
exploration . 

B.  APPLICATION  OF  THE  MERGING  THEOREM 

Next,  we  proceed  to  determine  how  well  the  merging 
theorem  predicts  the  timing  of  merging  events  was  performed. 
The  same  series  of  numerical  experiments  which  was  used  to 
verify  merging  of  layers  was  used  in  this  analysis.  As 
expected,  the  value  of  heat  flux  decreased  as  the  value  of 
Rp  increased  as  presented  in  Table  1.  Next,  the  magnitudes 
of  non-dimensional  heat  fluxes  FT  were  plotted  as  a  function 
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of  ,  and  the  curve  fitting  was  performed  using  Matlab 

curve  fitting  tool  in  order  to  determine  the  formula 
describing  the  coefficient  C  in  formula  (35)  .  Based  on  the 
best  fit,  a  formula  was  obtained 


C  =  6 . 2exp (-3.17  R^  ) 

Substituting  (36)  into  (35)  yields 

k  153exp(-6.24^)| 

B  9.5R  +4 

r  > 


(36) 


(37) 


This  form  of  XB  predicts  the  time  scale  for  a  particular 
merging  event  for  the  numerical  experiments  with  an  initial 
step  size  of  0.75.  The  inverse  of  XB  gives  time  in  non- 
dimensional  units.  The  time  scale  was  computed  using  formula 
(37)  and  compared  to  the  time  scale  determined  from 
numerical  experiments.  The  results  are  summarized  in  Table 
5. 


K 

theory 

4 

actual 

t 

theory 
(non-dim. ) 

t 

actual 
(non-dim . ) 

t 

actual 
(days ) 

1 . 1 

0 .0111 

0.01 

90.1 

100 

O 

00 

1.2 

0.0056 

0.008 

178.5 

125 

1 . 0 

1.3 

0.0028 

0.0029 

357 . 1 

340 

2.7 

1.4 

0.0014 

0.0011 

714 . 3 

910 

7 . 1 

Table  5.  Comparison  of  merging  theorem  with  numerical 

experiment . 
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In  analyzing  the  values  in  Table  5,  it  is  evident  that 
the  merging  theorem  predicts  the  time  scale  of  merging 
events  reasonably  well,  always  within  the  order  of 
magnitude.  The  verification  of  merging  theorem  for  higher 
values  of  Rp  was  not  performed  due  to  the  computational 

constraints;  each  increase  of  Rp  by  one  tenth  of  a  point 

results  in  an  approximate  doubling  of  the  time  scale.  Future 
studies  should  expand  the  verification  of  merging  theorem 
for  the  higher  values  of  R  . 


1.5 


1.4 


1.3 

Q. 

cc 

1.2 


1.1 


1 


30  60  90  100  130 

t 

Figure  21.  Variation  of  Rp  across  two  different  interfaces 

in  a  two-step  numerical  experiment.  The  green  (blue) 
curve  represents  the  higher  (lower)  initial  value  of 


The  physical  principles  described  in  Chapter  IV 


confirmed  the  results  in  Figure  21,  where  we  plot  the  same 
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records  of  the  density  ratios  at  two  interfaces  in  a  two- 
step  experiment  with  the  overall  R  =1 . 3  and  H=3  m.  Note  the 

continuous  increase  (decrease)  of  the  green  (blue)  curve 
representing  the  higher  (lower)  value  of  the  density  ratio 
at  t=60  after  the  initial  fluctuations.  Also  note  a  sharper 
decrease  in  the  lower  value  of  at  t=100  represented  by 

the  blue  curve  until  the  curve  reaches  R  =1,  the  moment 
when  merging  occurs . 
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V.  RESULTS  AND  CONCLUSIONS 


The  double-diffusive  process  has  received  a 
considerable  amount  of  attention  due  to  its  effect  on 
various  large  oceanic  features.  This  study  is  mainly 
focused  on  the  dynamics  of  diffusive  convection, 
particularly  the  numerical  modeling  of  the  dynamics  of 
thermohaline  staircases  and  associated  diffusive  fluxes. 
This  study  uses  the  results  of  the  study  of  oceanic  data 
from  the  Beaufort  Gyre  (Wilson,  2007)  in  order  to  compare 
the  magnitude  of  the  diffusive  heat  flux  and  to  analyze  its 
significance  on  the  Arctic  thermal  budget. 

In  order  to  compare  the  numerically  computed  diffusive 
heat  fluxes  to  the  fluxes  computed  from  oceanic  data,  the 
4/3  flux  is  invoked  (Turner,  1973) .  The  validity  of  the  4/3 
flux  was  tested  first  for  the  numerical  solution.  The 
results  are  summarized  in  Table  3.  The  analysis  revealed 
that  the  4/3  flux  law  works  reasonably  well  for  the  steps  < 
3  m,  but  the  amplitude  needs  adjustment. 

Padman  (1994)  and  Muench  et  al  (1990)  concluded  that 
diffusive  convection  could  play  a  significant  role  in 
influencing  the  vertical  heat  fluxes  from  the  warm  and  salty 
intermediate  water  upwards  to  the  polar  mixed  layers  and  sea 
ice.  However,  the  results  from  the  laboratory  experiments 
put  the  value  of  the  vertical  diffusive  heat  flux  on  the 
order  of  10~2  to  10_1,  a  non-significant  value  in  the  total 
Arctic  heat  budget.  Nevertheless,  it  is  necessary  to 
consider  that  the  laboratory  experiments  differ  in  many  ways 
from  the  actual  oceanic  conditions.  Diffusive  heat  fluxes 
computed  from  oceanic  data  from  the  Beaufort  Gyre  range  from 
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2  to  6  W  m2  (Wilson,  2007),  a  value  which  is  large  enough 
to  play  a  significant  role  in  the  Arctic  heat  budget. 

The  computed  values  of  the  diffusive  heat  fluxes  range 
from  1  to  3  W  rrT2  for  3<Rp<5  and  are  summarized  in  Table  4 

These  values  are  comparable  to  the  heat  fluxes  computed  from 
the  actual  oceanic  data.  However,  the  values  of  numerically 
computed  heat  fluxes  are  up  to  an  order  of  magnitude  greater 
than  the  results  from  the  laboratory  solutions.  Since  there 
is  large  discrepancy  among  these  results,  more  study  is 
necessary  before  the  final  assessment  of  the  role  of  the 
diffusive  vertical  fluxes  in  the  maintenance  of  the  Arctic 
thermocline  can  be  made. 

In  addition  to  computations  of  diffusive  fluxes,  this 
study  offers  some  insight  into  the  dynamics  of  thermohaline 
staircases.  The  merging  theorem  (Radko  2007)  predicts  types 
and  scales  of  merging  realized  in  the  thermohaline 
staircases.  The  results  of  this  study  confirm  that 
interfaces  within  the  diffusive  staircase  erode  over  time 
without  any  vertical  motion,  confirming  the  merging  theorem 
prediction  that  only  B-type  merger  is  taking  place. 

The  quantitative  predictions  of  the  merging  theorem 
were  tested  as  well  and  are  summarized  in  Table  5.  Based  on 
these  results,  the  numerical  experiments  confirm  that  the 
merging  theorem  predicts  the  time  scale  of  merging  events 
within  an  order  of  magnitude.  Due  to  the  computational 
constraints,  only  verification  for  the  lower  values  of  R^ 
was  performed.  Since  diffusive  convection  is  favorable  for 
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the  regions  with  2<R^<10  and  2<H<7  m,  future  studies  should 

focus  on  a  more  complete  exploration  of  this  parameter 
range . 
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